Thermodynamic properties and correlation functions of Ar films on the surface of a 

bundle of nanotubes 



O 

o 

(N 

u , 
Oh! 

< 



Nathan M. UrbanQ Silvina M. GaticaQ and Milton W. ColcH 

Department of Physics, Pennsylvania State University, 
University Park, Pennsylvania 16802-6300, USA 

Jose L. RiccardcH 

Departamento de Fisica, Universidad Nacional de San Luis, 5700 San Luis, Argentina 

We employ canonical Monte Carlo simulations to explore the properties of an Ar film adsorbed on 
the external surface of a bundle of carbon nanotubes. The study is concerned primarily with three 
properties: specific heat c(r), differential heat of adsorption qd, and Ar-Ar correlation functions 
g{r). These measurable functions exhibit information about the dependence of film structure on 
coverage and temperature. 
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I. INTRODUCTION 

Considerable interest has been attracted recently 
to the properties of simple gases (noble gases and 
small molecules) adsorbed near bundles of carbon 
nanotubesi^^r'^i4r^'^i7iVd"ii^^2,i'^i^4,i5,i6,i7,i8 ^his sub- 
ject has been reviewed recentlyii2*i^ The adsorption of 
these gases can occur within the tubes only if they are 
open, which is possible either during the process of nan- 
otube formation (e.g., when endohedral Ceo molecules 
are formed),^- or after chemical treatment to open the 
tubes. ^^'^^ The presence, or absence, of interstitial chan- 
nel (IC) molecules is an open question in the case of an 
idealized bundle of identical tubes; there seems to be no 
doubt, however, that such IC adsorption occurs in lab- 
oratory samples of polydisperse tubesiSS In contrast to 
variability in adsorption at these sites, the adsorption of 
gas on the external surface of the bundle is a ubiquitous 
phenomenon, in which the film coverage increases with 
the pressure (P) of the coexisting gas. In that exohe- 
dral environment, an adatom is strongly attracted to the 
groove region between two neighboring tubes; there, the 
film forms a quasi-one-dimensional phase. Further ad- 
sorption at low temperature (T) is predicted to manifest 
a so-called three-stripe phase of gas aligned parallel to 
the groovesi^ At higher gas coverage {N), there occurs 
a two-dimensional monolayer phase, qualitatively analo- 
gous to that found on the graphite surfaceiSSiS At even 
higher coverage, a multilayer film grows as P increases. 
There is an upper limit of total film coverage, set by the 
bundle's curvature ^2^*2^ this limit has yet to be explored. 

This study extends a previous investigation^ of the 
adsorption of Ar gas on the external surface of a nan- 
otube bundle. Argon was chosen as a model adsorbate 
because its gas-gas interaction is well known, making it 
a standard fluid in the study of simple fluids. In the 
previous study, denoted I, we employed the grand canon- 
ical Monte Carlo simulation method to explore the evo- 
lution of the equilibrium film as a function of P and T. 
The present paper, stimulated by recent and proposed 



experiments, adds three results to those derived in the 
previous study. One property is the specific heat, c(T), 
which is computed here from energy fluctuations, eval- 
uated using simulations within the canonical ensemble. 
The second property is the differential heat of adsorp- 
tion, Qd = —{dE/dN)T, where E is the energy of the 
film. This quantity is closely related to another quantity, 
which is more often measured experimentally, the isos- 
teric heat qst = {d{\nP) / dl3)N [where (3 = l/{kBT)], by 
the relationSi g^j = (/d + ksT (assuming an ideal gas co- 
existing with the film) . The third property reported here 
is the anisotropic correlation function of the overlayer. 
This quantity is related by Fourier transform to results 
of diffraction experiments. With the exception of the 
isosteric heat calculated by Shi and Johnson^SS none of 
these properties has been explored in simulation studies 
of films on nanotube bundles, prior to the present work. 

The outline of this paper is the following. Section UTI 
summarizes our simulation methods. Section IlIII reports 
results for the density and correlation functions. Sec- 
tion IIVI presents results for the thermodynamic proper- 
ties, c and qd- Section IVI summarizes our results. 



II. COMPUTATIONAL METHODS 

When not explicitly contradicted in this paper, it may 
be assumed that the physical system and computational 
method are as described in I. The primary model system 
is a bundle of infinitely long, cylindrically symmetric car- 
bon nanotubes of identical radii equal to 6.9 A. Only two 
adjacent nanotubes on the external surface of the bundle 
are simulated. The y axis is parallel to the nanotubes, 
and the z axis is directed away from the surface of the 
bundle. Periodic boundary conditions are imposed in the 
X and y directions (approximating the surface of the bun- 
dle as an infinite plane of nanotubes); reflecting bound- 
ary conditions arc imposed in the z direction. The unit 
simulation cell, whose volume contains half of each of the 
two adjacent nanotubes with the groove in between them 
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at the center, is 17 A in the x direction, 34 A in the y 
direction, and 40 A in the z direction. 

The simulations were done in the canonical ensem- 
ble, rather than the grand canonical ensemble more 
commonly used in adsorption simulations, in order to 
facilitate the calculation of the heat capacity. Two 
Markov Chain Monte Carlo simulation methods were 
used, the Metropolis algorithm.^^j^ and the Wang- 
Landau algorithmi^2i2ii2^ The Metropolis algorithm was 
used to calculate configurational observables, such as 
density distributions and correlation functions. The 
Wang-Landau algorithm was used to calculate thermo- 
dynamic observables expressible in terms of ensemble av- 
erages or their derivatives, such as specific and isosteric 
heat, for certain iV; the Metropolis algorithm was used 
to determine the full N dependence. 

The Metropolis algorithm proposes new configura- 
tions and accepts them with a probability equal to 
min {1, P(x')/P(x)}, where P(x) and P(x') are the prob- 
abilities of the old and new configurations x and x'; 
this acceptance rule causes the random walk to con- 
verge to the probability distribution -P(x). By choosing 
P(x) proportional to the Boltzmann factor exp [— /3[/(x)], 
where ?7(x) is the potential energy of the configuration, 
the Metropolis algorithm uniformly samples configura- 
tion space. For the Metropolis simulations of each (A^, T), 
4 X 10^ Monte Carlo moves were discarded during the ini- 
tial equilibration to converge the algorithm to the Boltz- 
mann distribution, then 4 x 10^ moves were generated, 
from which 10^ samples were drawn to perform simulated 
measurements of system observables. 

The Wang-Landau algorithm, like the Metropolis al- 
gorithm, also proposes and accepts configurations with 
a probability equal to min{l, P(x')/P(x)}. However, it 
chooses P(x) to be proportional to P[C/(x)] — l/(7[C/(x)], 
where g(U^ is the (relative) density of states, thus uni- 
formly sampling energy space (instead of configuration 
space, as in the Metropolis algorithm). It dynamically 
refines its estimate of the density of states by counting 
each visit to a state of a given energy U (or, rather, within 
a small range of energies U G [Ui — e/2,Ui + e/2] about 
an energy bin Ui of width e), and multiplies its running 
estimate of g{Ui) by a constant factor /. It continues the 
random walk until each energy is visited approximately 
uniformly (a "fiat histogram" of visits in energy space) , 
whereupon it reduces the factor / j^l"^ and starts 
another iteration. The algorithm terminates when / is 
reduced to a preset minimum greater than unity, with 
values closer to unity yielding more accurate estimates 
of the density of states. 

Once an estimate of g(C/) is produced, it can be used 
to calculate the partition function directly 

Z^|dxe-''^W«^g(C/.)e-^^^. (1) 

i 

Thermodynamic quantities can then be calculated from 
the partition function, as usual. One advantage of the 
Wang-Landau algorithm over the Metropolis algorithm 



(and the main reason for using it for this study) is that, 
because temperature dependence appears only in the 
Boltzmann weight exp(— /?[/) and not in the density of 
states (?([/) itself, a single simulation of g(t/) can calcu- 
late thermodynamic observables for all temperatures at 
once. 

Some modifications and improvements to the original 
published Wang-Landau algorithm were implemented. 
Boundary effects were properly handledi^^ To adapt the 
original lattice-based algorithm to continuum systems, 
preliminary Metropolis runs at low temperature were 
performed to estimate a lower bound on the energy bins 
(i.e., the ground state energy) The simulation can also 
become trapped for long periods of time in regions of 
high degeneracy, so that energies with small .g([/) go a 
long time before being revisited. To remedy this, the en- 
ergy bins can be broken up into overlapping subranges; 
ergodicity can be achieved more rapidly if the interval of 
energies to be traversed is smaller. Separate simulations 
are performed in each subrange, producing independent 
estimates of g{U). Some care must be taken in combining 
them into an estimate of g{U) over the full energy range: 
because each simulation calculates only the relative den- 
sity of states, the estimates will not generally match up 
at the boundaries of the subranges. To overcome this, 
each subrange estimate of g{U) is rescaled by a constant 
factor that minimizes the least-square error in log g(U) 
wherever two neighboring subranges overlap in energy^ 
This corresponds to choosing normalizing factors C„ that 
minimize the sum X^Jlog [9n{Ui)/Cn] - log [gn-i{Ui)]}'^ 
over the overlapping bins Ui (where gn denotes the den- 
sity of states simulated over subrange n) , and then rescal- 
ing gn{U) by C„. 

For the Wang-Landau simulations of each N, 1500 
equal-sized energy bins were used in a range [?7min,0], 
where ?7min is the ground state energy. The 1500 bins 
were divided into four overlapping subranges, simulated 
separately, consisting of the bins numbered 1-150, 76- 
787, 713-1425, and 1351-1500. A histogram was con- 
sidered "flat" when the number of visits to any partic- 
ular energy bin was less than ±20% the average num- 
ber of visits to any bin. The minimum / factor was 

/„,in = 1 + 10-^ 



III. CORRELATION FUNCTIONS 

For the purposes of this paper, the three-dimensional 
pair correlation function is defined as the probability den- 
sity g{r) that two particles are separated by a relative 
displacement r. Its projection G{x,y) = J dz g{r) into 
the xy plane is depicted in Fig. ^ The contours become 
wider and more irregular at higher temperatures, as the 
particles are thermally excited out of their well-defined 
low temperature sites. 

In the top pair of panels, one observes the nearly 
periodic, quasi-one-dimensional (ID) order within the 
groove. As studied recently in connection to nanotube 
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FIG. 1: Equiprobability contour plots of the projection into 
the xy plane of the pair correlation function, G{x, y) for, from 
top to bottom, = 9 (groove), N = I?,, N = 27 (three- 
stripe), and A'^ = 54 (monolayer), at T = 60 K (left) and 
T = 90 K (right). Distances are in angstroms. Note that the 
vertical axis scale is compressed relative to the horizontal. 



adsorption(2Si2S*2L224^2ii£ this phase may undergo a phase 
transition due to the weak interactions between particles 
in neighboring grooves. In the second pair of panels, 
one observes that the correlations within the three-stripe 
phase are weaker and even more T dependent than those 
in the groove phase. At 90 K, at higher coverage (seen in 
the middle four panels), the stripes are not as straight, 
primarily due to transverse excitation (as discussed in 
Sec. llVU|) . Note that the half-filled stripe case {N = 18) 
is somewhat less ordered than the completely filled three- 
stripe case {N = 27), as is expected. The bottom pan- 
els of Fig. n exhibit a highly correlated anisotropic two- 
dimensional (2D) solid at 60 K, the order of which washes 
out nearly completely by 90 K, as the monolayer melts. 

Previous experimental studies of Ar adsorption onto 
planar graphit6^ii^242ii4 found that the melting temper- 
ature depends on density, starting near 55 K at low den- 
sity, and increasing with density. Although the system 
studied here differs from that experiment in geometry, we 
expect the melting temperature of the monolayer in our 
system to similarly increase with density. The corruga- 
tion of the nanotube bundle should elevate the melting 
temperature somewhat compared to a planar surface, as 
the grooves will serve to more strongly confine the film's 
structure. 

Bienfait et al. have measured diffraction patterns for 



Ar on nanotube bundles.— Probably due to heterogeneity 
(for which there exists evidence in bare surface diffrac- 
tion), the diffraction data are not easily interpreted. 
There is, however, definite evidence of ID interatomic 
spacing (i.e., a peak at wave vector q = 17/nm) at 
low coverage and 2D close-packed spacing (peak near 
20/nm). 



IV. HEAT CAPACITY 

A. Overview 

The isochoric specific heat, c{T) = {dE/dT)v/N, 
i.e., the heat capacity (per particle) at constant volume, 
was calculated from ensemble averages. It is known^ 
that the heat capacity can be given in terms of en- 
ergy fluctuations, {dE/dT)v = ((-B^) - {Ef)/{kBT^), 
where (•) denotes an expectation taken over the canon- 
ical ensemble. The equipartition theorem gives the ki- 
netic energy contribution of ^ks per degree of free- 
dom to the specific heat, yielding a total c(T) = 

ffcs + (([/2) - {uf)/{NkBT^). Given the density 
of states g{U) calculated with the Wang-Landau algo- 
rithm, the expectations may be calculated from ([/) — 
Y.i Ui g{Ei) exp (-/3C/,), and similarly for (C/^). The 
heat capacity was also estimated directly from the deriva- 
tive {dE/dT)v by means of a finite difference approxi- 
mation. These latter estimates, while consistent with the 
fiuctuation estimates, were "noisier" and are not consid- 
ered further in this paper. 

The simulated c(T) for the groove, three-stripe, and 
monolayer phases is shown in Fig. |2| Note that the over- 
all trend is for c{T) to have a remarkably high value, 
in the range 3-7 Boltzmanns, much higher than might 
be expected from simple quasi-one-dimensional and two- 
dimensional models. We do not have a detailed quantita- 
tive model to explain all of the observed features, but we 
can give a qualitative explanation of its behavior. The 
explanations are justified by examining the probabilities 
of finding particles at given energies in the external po- 
tential. Fig. 13 indicating the fraction of particles that 
are in the groove, stripes/monolayer, etc. (quantified in 
Table P). 



B. Low density 

Consider first the low-density limit. At low temper- 
atures, the specific heat is near 2.5 Boltzmanns. This is 
to be expected: the three kinetic degrees of freedom each 
contribute the usual 1/2 Boltzmann; the two transverse 
dimensions, for which the external potential is approx- 
imately harmonic at its minimum at the center of the 
groove, each contribute another 1/2 Boltzmann. As the 
temperature increases, a peak in the specific heat oc- 
curs near 170 K when substantial numbers of adatoms 
are promoted out of the groove and into monolayer sites 
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FIG. 2: The dimensionless specific heat c{T)/kB of the low- 
density limit (A'' = 1); the groove (A'^ = 9), three-stripe (A^ = 
27), and monolayer (A'^ = 54) phases; and the theoretical 
prediction of the low-density limit given by the dimensional 
crossover model (discussed in the text). 
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FIG. 3: The relative probability P(t/cxt) of finding a parti- 
cle at an external potential of L^ext, for the low-density limit 
[N = 1), and the groove (A^ = 9), three-stripe (A^ = 27), and 
monolayer (A'' = 54) phases, for various temperatures. 
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TABLE I: Percentage of particles in the groove ( — 1600 K < 
(7oxt < -1200 K), monolayer (-1200 K < C/cxt < -400 K), 
and vapor ((7ext > —400 K) regimes, for the low-density limit 
(A' = 1), and the groove (A^ = 9), three-stripe (A' = 27), and 
monolayer (A' = 54) phases, for various temperatures. 



limit, can also be understood by examining the so- 
called "volume density of states"^ f{U), defined such 
that /([/) dU is the volume of space bounded by in- 
finitesimally separated isopotential surfaces, U < Ucxt < 
U + dU. This function is related to the energy proba- 
bility density P{U) in Fig. O at low density by P{U) — 
pf{U) exp {—(3U), where p is the number density of parti- 
cles. By dividing an estimate of P{U) at a given temper- 
ature (here, T = 300 K) by the exponential Boltzmann 
factor, we obtain an estimate proportional to the volume 
density of states /(C/), depicted in Fig. ^ 




-1600 -1400 -1200 -1000 -800 -600 -400 -200 
U (external) [K] 



elsewhere on the surface of the nanotubes (see Table Pi . 
As r ^ oo, the adatoms desorb from the surface alto- 
gether, and c(T) approaches the 3/2 Boltzmanns of the 
three kinetic degrees of freedom of a pure vapor. (This 
will be the case for all other densities as well, in the high- 
T limit.) 

These conclusions are corroborated, as mentioned, in 
the first row of Fig. O these results, in the low-density 



FIG. 4: The function f{U) (unsealed), giving the volume of 
space f{U) dU enclosed within a range of external potential 
energy [U, U + dU]. 

The qualitative form of this figure can be explained 
by appealing to a previously studied analytic model, 
the dimensional crossover modeL^ This exactly soluble 
model ignores interparticle interactions (an assumption 
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appropriate for the low-density limit) and treats the nan- 
otube bundle as consisting solely of two regions, a one- 
dimensional groove region approximated by a harmonic 
potential in the two transverse dimensions, and a two- 
dimensional planar monolayer region approximated by a 
harmonic potential normal to the surface. Evaporation 
from the monolayer to vapor is neglected. Its (configura- 
tional) partition function is given by 

- f w2„p-/3[Vg-H(l/2)ar=] 
■^crossover / u- / c- 

J groove 

+ lJ dze-''[^'"+(i/2)'='"-'l, (2) 

<y mono 

where the parameters Vg = -1671 K, a = 4898 K/A^, 
Vm = -853 K, and fc™ = 4792 K/A^ were determined by 
a fit to the external potential, = 18 A is the approx- 
imate width of the monolayer region in the transverse 
direction, and integrations were taken over regions ex- 
tending 2 A away from the groove minimum and 1 A 
away from the monolayer minimum. 

For the low energies dominated by the groove phase, 
the crossover model approximates the external potential 
as J7 = -I- ^ar^. The cylindrical volume enclosed by 
an isopotential goes like V ^ r'^, and f{U) = dV/dU = 
{dV/dr)/{dU/dr), which is a constant; indeed, the f{U) 
calculated in Fig. 21 is nearly constant at low energies. 
For the monolayer region, close to the potential mini- 
mum, U — Vm + ^kmZ^ ■ The rectangular volume en- 
closed by an isopotential goes like V ^ z, and f{U) — 
{dV/dz)/{dU/dz), which goes like z-^ [U - VmY^''^ 
for [/ > Vm- This divergence in the model accounts 
qualitatively for the peak in /(J7) just above the mono- 
layer energy of about —800 K. For high energies dom- 
inated by the vapor phase, we can treat the substrate 
as a semi-infinite rectangular volume, and approximate 
the external potential by a long-distance Lennard- 
Jones potential integrated over this region, which yields 

U Then /([/) wiU go like z"* - {-Uy^l^ , 

qualitatively accounting for the sharp rise in /([/) as 
C/^0. 

Given the volume density of states f{U), we can then 
obtain the energy probability density P{U) at any other 
temperature by scaling this temperature-independent 
function by the appropriate Boltzmann weight. In par- 
ticular, the three columns of the first row of Fig. O are 
just the function in Fig. 0] scaled by Boltzmann factors 
exp {—f3U) that decay with decreasing rapidity as T in- 
creases {P decreases). The groove is highly populated at 
low temperatures (large /?) when the exponential damp- 
ing is great enough to suppress population at higher ener- 
gies. The monolayer becomes populated at intermediate 
temperatures when the damping is no longer sufficient 
to suppress the peak in f{U) that occurs at the mono- 
layer energy (—800 K), and the vapor becomes populated 
at still higher temperatures (small /3) when the damping 
fails to suppress the rapid increase in f{U) towards K. 

The partition function of the crossover model can also 



be used to calculate the specific heat directly. As seen 
in Fig. [3 the correspondence between the prediction of 
this analytic approximate model and the simulated full 
model is quite good. 



C. Higher coverage 

Next, consider the groove phase. At low temperatures, 
the specific heat is near 3 Boltzmanns. A contribution 
of 2.5 Boltzmanns is accounted for by the same argu- 
ment as for the low density limit. Unlike the low density 
limit, however, the groove phase is densely packed with 
adatoms, and interparticle interactions must be consid- 
ered. An additional 1/2 Boltzmann arises from confine- 
ment in the longitudinal dimension, for which the Ar- 
Ar interaction potential is approximately harmonic at 
its minimum when the adatoms are stably distributed in 
equilibrium. As the temperature increases, evaporation 
out of the groove begins. At T ss 70 K, evaporation is 
great enough to excite adatoms out of the groove; while 
not many of these atoms reach the monolayer region (Ta- 
ble P), there is still a large change in potential energy for 
a small increase in temperature, and thus a large specific 
heat. The specific heat then decreases slightly with in- 
creasing temperature, since the change in energy is not 
as large once the initial adatoms have begun to be pro- 
moted. This low-temperature peak is not present in the 
low-density case because, as seen in both Table and 
Fig.n the adatoms are not spread transversely as greatly 
about the immediate groove region at low temperatures 
in the low-density case as they are in the groove case. 
However, in a manner qualitatively analogous to the low- 
density limit, an additional, larger peak in the specific 
heat is found at still higher temperatures (T w 175 K), 
mostly from promotion from the groove into the stripes 
and the rest of the monolayer. 

Like the groove phase, the three-stripe phase starts 
out at 3 Boltzmanns at low temperatures, similar to the 
groove phase, except that the 1 Boltzmann from external 
potential confinement in the transverse plane is replaced 
by 1/2 Boltzmann from external potential confinement 
normal to the surface, and 1/2 Boltzmann from inter- 
particle confinement along the surface in the transverse 
plane. The specific heat peaks near T = 55 K; this is 
not due to a significant fraction of particles being pro- 
moted from the groove to the stripes, as one might ex- 
pect, but rather to a wider range of energies within the 
stripe/monolayer region, and promotion from the stripes 
to the rest of the monolayer; see Fig. Oand Table It 
peaks again near T — 165 K, with as the groove empties 
into the stripes and monolayer, as well as the beginning 
of evaporation off the surface. 

The monolayer phase also starts out at 3 Boltzmanns 
at low temperatures, for reasons analogous to the three- 
stripe phase. At higher temperatures, there is a peak in 
the specific heat near T = 80 K which, like the three- 
stripe peak, is largely due to a broadening of the parti- 
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cles across a range of energies in the monolayer region, 
as well as some promotion from the monolayer to the bi- 
layer (Fig. EJ. Another peak appears near T — 140 K, 
corresponding to evaporation out of the monolayer into 
the bilayer, and to vapor. 



D. Further results 

The N dependence of several isotherms is displayed in 
Fig. Of particular note is the rapid rise in the specific 
heat near TV ~ 8, just before groove completion, as T 
goes from 60 to 90 K. This is attributed to promotion 
out of the groove. Similarly, near monolayer completion 
the marked increase in c with T is attributed to thermal 
promotion out of the monolayer. 




10 15 20 25 30 35 40 45 50 55 
N 



FIG. 5: The dimensionless specific heat c{T) /fcs as a function 
of density, at T = 60, 75, and 90 K. 

It is also illuminating to study the differential heat of 
adsorption, qd{N) — —{dE/dN)T, the energy required 
to adsorb an additional particle onto the surface at con- 
stant temperature. The differential heat is related to 
the heat capacity at constant density, Cn — {dE/dT)^, 
by a Maxwell relation derived from the total derivative 
dE = {dE/dT)N dT + [dE/dN)T dN, which yields 



dC, 



N 



dN 



dT 



(3) 



N 



(Note that Cn — Cy in the canonical ensemble.) 

The differential heat of adsorption is summarized in 
Fig. El At low densities, the differential heat is near the 
minimum of the external groove potential « —1600 K, 
becoming slightly larger at lower temperatures. Both 
this value and the T dependence at low N can be un- 
derstood from the low-density equation of state, U/N — 



^ksT. As additional particles are added, the en- 



As the groove phase is approached, the groove becomes 
tightly packed and the interaction energy becomes signif- 
icant, so that adding an additional particle reduces the 
energy by the external groove potential plus the Lennard- 
Jones well depth, e « —120 K for Ar. The difference in 
QdiN) between low and high temperatures is particularly 
large just before the groove phase, which in accordance 
with Eq. (O corresponds to the steepest increase in heat 
capacity as seen in Fig. [SJ at the = 9 groove phase 
itself, where the heat capacity peaks with increasing N, 
we see correspondingly little difference in the differen- 
tial heat at various temperatures. At the other extreme, 
near monolayer completion, a similar T dependence is 
observed. The large decrease in qd with increasing T is 
consistent with Eq. © and Fig.|Sl the latter shows a large 
value of dC/dN except below 60 K. The explanation is 
monolayer-to-bilayer promotion above 60 K. 




G^T = 60K 
QQT = 75 K 
AAT = 90K 



ergy for each additional particle is reduced by slightly 
more than this amount, to include the interaction energy. 
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FIG. 6: The differential heat of adsorption qd{N)/kB at T — 
60, 75, and 90 K. [Shift each curve upwards by its temperature 
T to obtain the isosteric heat qst{N)/kB.] 

Experimental measurements of isosteric heat for ar- 
gon on nanotube bundles have been reported by Wil- 
son et al.^ Talapatra, Rawat, and MigonCfi Jakubek 
and Simard,"* and Bienfait et al.^ grand canonical 
Monte Carlo simulations have been published by Shi and 
Johnsoni^ 

The isosteric heat calculations of Shi and Johnson for 
adsorption of Ar on a homogeneous bundle at 90 K agree 
closely with our results, with a peak of qst = 14 kJ/mol 
just before the groove phase, corresponding to our peak of 
1650 K. Past the groove phase, their calculated isosteric 
heat drops and remains constant with coverage, slightly 
below 10 kJ/mol, corresponding to our nearly constant 
value near 1200 K. 

Shi and Johnson compared their calculations to the ex- 
perimental results of Wilson et al. and Talapatra et al. , 
and since we agree with those calculations, we will briefly 
summarize their conclusions. Our calculations agree with 
both experiments at higher coverage, beginning at the 
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three-stripe phase, but their isosteric heats at lower cov- 
erage are dramatically greater than ours, as large as 
18 kJ/mol ('^ 2200 K) at low coverage. We ascribe this 
discrepancy with experiment to our neglect of bundle het- 
erogeneity, following Shi and Johnson, whose simulations 
of heterogeneous bundles agreed well with both experi- 
ments. 

In contrast, the isosteric heats measured by Jakubek 
and Simard agreed well with our simulations, with a peak 
of 137 meV (~ 1600 K) near the groove, descending to 
plateau of 106 meV (~ 1200 K) through to monolayer 
coverage. This agreement with our calculations suggests 
that their bundles were more homogeneous than those 
studied in the other two experiments. It should be noted, 
however, that their isosteric heat continues to decrease as 
coverage increases, whereas our isosteric heat appears to 
rise slightly as the monolayer is approached. The isosteric 
heat of Wilson et al. also drops past the monolayer. 

Like the other experiments, the results of Bienfait et al. 
for Ar exhibit two plateaus in the dependence of the isos- 
teric heat on coverage. The lowest coverage data yield 
qst = 15 kJ/mol, or about 1800 K. Our predicted value 
in this range is of order 1650 K. The higher coverage, 
broad plateau corresponds to a measured qst — 1200 K, 
which agrees well with the value we find for the three- 
stripe phase. However, the data at monolayer coverage 
continue to decrease, while ours appear to increase, as 
noted. Another area of disagreement is the extent, in 
coverage, of these plateaus. In the data, the second 
plateau extends over a coverage range comparable to that 
of the first plateau. Our calculations, instead, find that 
groove region of high q^t extends over just one-sixth of 
the range of the combined three-stripe plus monolayer 
regime (grouped together because of similar values of 
qst)- This discrepancy may be attributed to the role of 
large interstitial cavities within the bundle, as argued by 
Bienfait et al. 



V. SUMMARY AND CONCLUSIONS 

Our results are intended to stimulate further exper- 
imental studies of this system and analogous systems 
involving other gases on nanotube bundles. We have 
investigated the variation of thermodynamic properties 
with T and N . One of the more interesting general 
results is that the specific heat is typically larger than 
might have been expected from either simple models 
used to treat these systems (either independent parti- 
cles or a solid) or from experimental results for films 
on graphite^ For most conditions studied here, the spe- 
cific heat exceeds three Boltzmanns, with average val- 
ues in the range four to five Boltzmanns. In contrast, 
the specific heat of independent particles^ in this envi- 
ronment is less than three Boltzmanns, except at high 
T (above 100 K), when the particles are excited out 
of the groove. The large values found in these simula- 
tions arise from the fact that the highly corrugated po- 



tential surface presents a sequence of excitation steps 
(groove— i-three-stripe—»monolayer—> • • ■ vapor), each 
of which enhances the specific heat. 

The temperature dependence of the specific heat shows 
a characteristic double-peak structure. All densities show 
a large peak near 175 K, corresponding to promotion of 
adatoms out of the groove into the monolayer region. 
The groove, three-stripe, and monolayer phases show an 
additional peak at lower temperature, corresponding to 
a thermal broadening in the range of external potential 
energies of the particles, rather than to any significant 
promotion of particles into qualitatively different regions. 

Other principal results involve the relation between the 
evolution of film structure (with increasing TV) and the 
corresponding thermodynamic and correlation functions. 
As the groove begins to fill {N approaching 9), the heat 
capacity shows a dramatic jump as a function of cover- 
age. Consistent with the Maxwell relation, Eq. the 
differential heat decreases with T at that point (Fig.|HJl. 
Analogous behavior occurs near completion of the three- 
stripe phase, near N = 27. 

Our study has been fully classical, but the tempera- 
tures beneath which quantum effects become significant 
can be estimated]^ We estimate that quantum effects 
can be ignored above about 80 K (see Appendix 
This is a higher temperature than some of the impor- 
tant structure in the heat capacity — the first peak in the 
heat capacity occurs at or below this temperature. Mod- 
ifications to the heat capacity from quantum mechanics 
at very low energies are given by Debye theory 
we expect that c(T) ^ as T ^ 0, and c(T) cx T"^, where 
d ~ 1 for the groove and d ^ 2 for the monolayer, if the 
density is high enough to form a bulk phase. To evaluate 
quantum effects accurately would require application of 
the path integral Monte Carlo method to the problem. 

We note, also, that an experimental heat capacity cell 
has a volume on the order of 1-10 cm^, whereas our sim- 
ulation volume was on the order of lO"^'' cm^. Our sim- 
ulation, focusing on small volume nearer the adsorbed 
film, thus ignores almost all of the volume in which des- 
orption into vapor can occur. This causes the simulation 
to underestimate the heat capacity that will be experi- 
mentally measured. The effects of desorption cannot be 
ignored when the number of atoms in the vapor starts to 
approach the number of atoms in the film; this occurs at 
roug hly 25-50 K (see Appendix lU . 

Particularly interesting results from the correlation 
function studies include the reduced longitudinal corre- 
lations in the groove and striped phases as T rises above 
60 K. These results would be amenable to testing by 
diffraction experiments even if the samples included a 
randomly oriented batch of nanotubes; this is a familiar 
problem dealt with in powder averaging of small-sample 
experiments. 

This paper studied a system of identical nanotubes. 
The sensitivity of c(T) to nanotube heterogeneity, with 
an asymmetric groove region between nanotubes of dif- 
ferent sizes, is a potentially interesting subject for future 
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nificant, by determining when the ratio Ny/Nm of atoms 
in the vapor to atoms in the monolayer becomes signif- 
icant. Extending the crossover model, we can consider 
the monolayer and vapor phases as separate systems, the 
monolayer modeled as a surface with a harmonic normal 
potential, and the vapor modeled as a free gas. The ra- 
tio Ny/Nm is then given by the ratio of their respective 
partition functions 



dz/ [ dze-'^[^'"+(i/2)fe™^^]. (Bl) 

m J cell ' J mono 



APPENDIX A: QUANTUM EFFECTS 

The upper bound on the temperature at which quan- 
tum effects must be considered is dominated by the 
physics of the deepest energy well, i.e., the groove. We 
can obtain one estimate by considering the minimum 
energy of longitudinal phonons in the groove in Debye 
theory, huo, where ujo — \fkfm and k = 2^/'^(9e/cr^) 
is the force constant of a quadratic approximation to 
the minimum of the Ar-Ar interaction potential, a 12- 
6 Lennard- Jones potential C/int = 4e[(CT/r)^^ — (a/r)^], 
with a = 3.4 A, e = 120 K for Ar. The correspond- 
ing energy is 27 K. This estimate considers only Ar-Ar 
interactions and ignores the external potential; we can 
obtain a complementary estimate by ignoring the inter- 
actions and considering only the external potential. We 
return to the crossover model of the groove, outlined in 
Sec. IIVI as a two-dimensional harmonic oscillator with 
a force constant a = 4898 K/A'^. Treating it now as a 
quantum harmonic oscillator, it is excited at an energy 
fiiLj±, where uj± = ^ a/m and m is the atomic mass of 
argon. The corresponding energy of this second estimate 
is 77 K. Taking the larger of the two as a conservative 
estimate, we expect that quantum effects can be ignored 
above about 80 K. 



APPENDIX B: EFFECTS OF DESORPTION 

We can estimate the temperature at which desorption 
into the full volume of an experimental cell becomes sig- 



We take the first integral between zero and the cell height, 
h] the second integral may be safely taken between zero 
and oo, in the interests of finding an analytic solution. 
This gives 



m } 



The fraction of atoms in the vapor above which desorp- 
tion should be considered "significant" is ambiguous, but 
we might take it to be 10%-20%. Solving Eq. for /3, 
using h = 1 cm and the values for Vm and km found in 
Sec ■ II V 51 this corresponds to a temperature in the range 
of 25-50 K. 

This calculation neglects interparticle interactions. 
Their inclusion would lower the estimate of the temper- 
ature at which desorption from the monolayer into va- 
por becomes significant, analogously to how the evapo- 
ration from the groove to the monolayer takes place at a 
lower temperature when the groove is packed — the inter- 
acting case — than when it is sparsely populated and the 
adatoms are effectively noninteracting. This is supported 
by the data in Table ^ rnore groove^monolayer promo- 
tion occurs for A^ — 9atr=175K than for iV = 1 at 
the comparable temperature T = 170 K, indicating that 
the groove promotion begins at a lower temperatures for 
the interacting = 9 than for the noninteracting = 1 . 
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